#!/bin/csh
#...........................................................................
#BSUB -a poe                    # use LSF poe elim
##BSUB -x                        # exclusive use of node (not_shared)
#BSUB -n 1                      # total tasks needed
#BSUB -R "span[ptile=16]"       # max number of tasks (MPI) per node
#BSUB -J NDGV1
#BSUB -o %J.runout              # output filename
#BSUB -e %J.runout              # error filename
#BSUB -W 05:59
##BSUB -P P93300642
#BSUB -P WYOM0024
##BSUB -q regular             # queue
##BSUB -q small                  # regular 
#BSUB -q caldera

set echo 

# CAM5 orography 
# 
# /glade/p/cesm/cseg/inputdata/atm/cam/topo/USGS-gtopo30_1.9x2.5_remap_c050602.nc oro_ecmwf_T170.nc
#  
# cdo sp2gp -sp2sp,170 -gp2sp -remapcon,t170grid oro_cam5.nc oro_cam5_t170.nc


foreach yy (2013) 
###foreach mm (01 02 03 04 05 06 07 08 09 10 11 12) 
foreach mm (10) 
foreach dd (01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31) 

setenv mpath /glade/p/rda/data/ds627.0/ei.oper.an.ml/${yy}${mm}
setenv tpath /glade/scratch/kaizhang/data/ndata/tmpv1
setenv opath /glade/scratch/kaizhang/data/ndata/output
setenv ipath /glade/scratch/kaizhang/data/ndata/input
setenv spath /glade/scratch/kaizhang/data/ndata/cam5nud

mkdir -p ${tpath} 

foreach hh (00 06 12 18) 

echo ''
echo ''
echo ${yy}${mm}${dd}${hh} 
echo ''
echo ''

rm -f ${tpath}/ei.oper.an.ml.regn128sc.${yy}${mm}${dd}${hh}.grb
rm -f ${tpath}/ei.oper.an.ml.regn128uv.${yy}${mm}${dd}${hh}.grb

ln -s ${mpath}/ei.oper.an.ml.regn128sc.${yy}${mm}${dd}${hh} ${tpath}/ei.oper.an.ml.regn128sc.${yy}${mm}${dd}${hh}.grb
ln -s ${mpath}/ei.oper.an.ml.regn128uv.${yy}${mm}${dd}${hh} ${tpath}/ei.oper.an.ml.regn128uv.${yy}${mm}${dd}${hh}.grb

#.............................................................................. 
# 3. Do vertical interpolation using cdo remapeta 
#.............................................................................. 

#       pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybi_cam(k)) + 1013.25 * (hyai_cam(k))
#
#       pold(k,:,:) = ( surp(itim,:,:)/100.) * (hybi_echam(k)) + (hyai_echam(k)) / 100  
#
 rm -f ${tpath}/tmp${yy}${mm}${dd}${hh}.grb 

 cdo merge ${tpath}/ei.oper.an.ml.regn128sc.${yy}${mm}${dd}${hh}.grb ${tpath}/ei.oper.an.ml.regn128uv.${yy}${mm}${dd}${hh}.grb ${tpath}/tmp${yy}${mm}${dd}${hh}.grb 

 rm -f ${tpath}/tmp${yy}${mm}${dd}${hh}_stps.nc 
 rm -f ${tpath}/tmp${yy}${mm}${dd}${hh}_uuvv.nc 

 cdo remapeta,vct_L30,oro_cam5_t170.nc  ${tpath}/tmp${yy}${mm}${dd}${hh}.grb ${tpath}/tmp${yy}${mm}${dd}${hh}_stps.nc 
 cdo remapeta,vct_L30,oro_cam5_t170.nc  ${tpath}/tmp${yy}${mm}${dd}${hh}.grb ${tpath}/tmp${yy}${mm}${dd}${hh}_uuvv.nc

 rm -f ${tpath}/tmp${yy}${mm}${dd}${hh}_o1.nc 
 rm -f ${tpath}/tmp${yy}${mm}${dd}${hh}_o2.nc 

 cdo remapcon,grid_2x2 -expr,'var152=exp(var152);var130=var130;var156=var129/9.8;' -selvar,var130,var152,var129 ${tpath}/tmp${yy}${mm}${dd}${hh}_stps.nc ${tpath}/tmp${yy}${mm}${dd}${hh}_o1.nc 
 cdo remapcon,grid_2x2 -selvar,var131,var132 ${tpath}/tmp${yy}${mm}${dd}${hh}_uuvv.nc ${tpath}/tmp${yy}${mm}${dd}${hh}_o2.nc 

 rm -f ${opath}/ndg_v1_2x2.5_${yy}${mm}${dd}${hh}.nc 

 cdo merge ${tpath}/tmp${yy}${mm}${dd}${hh}_o1.nc ${tpath}/tmp${yy}${mm}${dd}${hh}_o2.nc ${opath}/ndg_v1_2x2.5_${yy}${mm}${dd}${hh}.nc 

#.............................................................................. 
# 3. Do vertical interpolation using cdo remapeta 
#.............................................................................. 

cat > make_nudge_z3_${yy}${mm}.ncl << EOF

  load "/glade/apps/opt/ncl/6.1.2/gnu/4.4.6/lib/ncarg/nclscripts/csm/contributed.ncl"

  fna = "${opath}/ndg_v1_2x2.5_${yy}${mm}${dd}${hh}.nc" 
  fno = "${opath}/z3_2x2.5_${yy}${mm}${dd}${hh}.nc"

  fla = addfile(fna,"r") 

  st = fla->var130 
  ps = fla->var152 

  ;;printVarSummary(st) 
  ;;printVarSummary(ps) 

  hyam = fla->hyam 
  hybm = fla->hybm 
  hyai = fla->hyai 
  hybi = fla->hybi 

  pmid = st 

  nk = 30 

  nlat = dimsizes(st(0,0,:,0)) 
  nlon = dimsizes(st(0,0,0,:)) 

  pint = new((/ 1, nk+1, nlat, nlon /), "float")  

  pint!0 = "time" 
  pint!1 = "lev" 
  pint!2 = "lat" 
  pint!3 = "lon" 

  ;;printVarSummary(pint) 

  ps0 = 1.e5 

  do ik = 0,nk-1 
     ;;;pmid(0,ik,:,:) = (/ dble2flt(hyam(ik)) + dble2flt(hybm(ik)) * exp(ps(0,:,:)) /)  
     pmid(0,ik,:,:) = (/ dble2flt(hyam(ik)) + dble2flt(hybm(ik)) * (ps(0,:,:)) /)  
  end do 

  print(" ") 
  print(" ") 
  print(" ") 
  print("---------- pmid --------------") 
  print(pmid(0,:,{0.},{0.}))  

  do ik = 0,nk
     ;;pint(0,ik,:,:) = (/ dble2flt(hyai(ik)) + dble2flt(hybi(ik)) * exp(ps(0,:,:)) /)    
     pint(0,ik,:,:) = (/ dble2flt(hyai(ik)) + dble2flt(hybi(ik)) * (ps(0,:,:)) /)    
  end do 


  pdel = pmid 

  do ik = 0,nk-1
     pdel(0,ik,:,:) = (/ pint(0,ik+1,:,:) - pint(0,ik,:,:) /) 
  end do

  print(" ") 
  print(" ") 
  print(" ") 
  print("---------- pdel --------------") 
  print(pdel(0,:,{0.},{0.}))  

  gravit = 9.8 
  rair = 6.02214 * 1.38065 * 1.e3 / 28.966
  ;;rair = 287.058 

  rog = rair/gravit

  zi = pint 
  zm = pmid 

  ;; The surface height is zero by definition.

  zi(0,nk,:,:) = 0.0

  
  hkk = st(0,0,:,:) 
  hkl = st(0,0,:,:) 

  ;;printVarSummary(pmid) 
  ;;printVarSummary(pdel) 
  
  ;; Compute zi, zm from bottom up. 
  ;; Note, zi(i,k) is the interface above zm(i,k)

  do k = 0, nk-1 

     ;; compute zm, zi

     kk = 29 - k 

     ;;print(" kk = : "+kk) 

     hkl(:,:) = (/ pdel(0,kk,:,:) / pmid(0,kk,:,:) /) 

     hkk(:,:) = (/ 0.5 * hkl(:,:) /) 

     zm(0,kk,:,:) = (/ zi(0,kk+1,:,:) + rog * st(0,kk,:,:) * hkk(:,:) /) 
     zi(0,kk,:,:) = (/ zi(0,kk+1,:,:) + rog * st(0,kk,:,:) * hkl(:,:) /) 

     ;;print( pdel(0,kk,:,:)+ "   "+pint(0,kk-1,:,:)+"   "+pint(0,kk,:,:)+"    "+zm(0,kk,:,:) +"    "+ zi(0,kk,:,:)+ "    "+ st(0,kk,:,:) + "   " +  hkk(:,:)) 

  end do

  ;;print(zm(0,29,0,20) + "    "+ zm(0,10,0,20) +"    "+  zm(0,1,0,20)) 


  system("rm "+fno) 

  flo = addfile(fno,"c") 

  flo->Z3=zm

EOF


 ncl make_nudge_z3_${yy}${mm}.ncl



#
# cdo remapcon,grid_2x2 ${tpath}/tmp${yy}${mm}${dd}${hh}_z3.nc ${opath}/z3_2x2.5_${yy}${mm}${dd}${hh}.nc
#

end # end hh 

end # end dd 

cdo mergetime ${opath}/ndg_v1_2x2.5_${yy}${mm}????.nc ${opath}/ndg_v1_2x2.5_${yy}${mm}.nc
mv  ${opath}/ndg_v1_2x2.5_${yy}${mm}????.nc bak 

cdo mergetime ${opath}/z3_2x2.5_${yy}${mm}????.nc ${opath}/z3_2x2.5_${yy}${mm}.nc
mv  ${opath}/z3_2x2.5_${yy}${mm}????.nc bak


end # end mm

end # end yy




